# ------------------------------------------------------------------------------#
# Description: Create Table 7: Peer effects heterogeneity by visibility and sign placement
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(pacman)
p_load(data.table, fixest, stargazer)
options(scipen = 999)

# Load and merge base data -----------------------------------------------------
load("data/adoptions/adoptions.RData")
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale     := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale  := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale      := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale      := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale    := peer.adopt.both.100 / 100]

dat <- dat[year > 2010]
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Load and merge GSV visibility and sign data ----------------------------------
load('data/gsv/rw_gsv.RData')
load('data/gsv/peer_adopt_buckets_any_vis.RData')
dat <- merge(dat, peer.adopt.bucket.any.vis, by = c("pin", "year"))
rm(peer.adopt.bucket.any.vis)

load('data/gsv/peer_adopt_buckets_any_rg_vis.RData')
dat <- merge(dat, peer.adopt.bucket.any.rg.vis, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg.vis)

load('data/gsv/peer_adopt_buckets_cs_vis.RData')
dat <- merge(dat, peer.adopt.bucket.cs.vis, by = c("pin", "year"))
rm(peer.adopt.bucket.cs.vis)

load('data/gsv/peer_adopt_buckets_any_sign.RData')
dat <- merge(dat, peer.adopt.bucket.any.sign, by = c("pin", "year"))
rm(peer.adopt.bucket.any.sign)

load('data/gsv/peer_adopt_buckets_any_rg_sign.RData')
dat <- merge(dat, peer.adopt.bucket.any.rg.sign, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg.sign)

load('data/gsv/peer_adopt_buckets_cs_sign.RData')
dat <- merge(dat, peer.adopt.bucket.cs.sign, by = c("pin", "year"))
rm(peer.adopt.bucket.cs.sign)

# Estimate models --------------------------------------------------------------
ivsign <- feols(adopt.any*100 ~ peers.100 | blkgrp + year |
                  peer.adopt.any.100 + peer.any.sign.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
                cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

ivvis <- feols(adopt.any*100 ~ peers.100 | blkgrp + year |
                 peer.adopt.any.100 + peer.any.vis.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
               cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

ivonlysign <- feols(adopt.any*100 ~ peers.100 | blkgrp + year |
                      peer.any.sign.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
                    cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

ivonlyvis <- feols(adopt.any*100 ~ peers.100 | blkgrp + year |
                     peer.any.vis.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
                   cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

ivboth <- feols(adopt.any*100 ~ peers.100 | blkgrp + year |
                  peer.any.sign.100 + peer.any.vis.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
                cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

# Define variable labels -------------------------------------------------------
varnames <- c(
  'year' = 'Year',
  'blkgrp' = 'Block Group',
  'elig.peers.cs.100' = '\\# Eligible Peers (CS)',
  'elig.peers.rg.100' = '\\# Eligible Peers (RG)',
  'peer.adopt.any.100' = '$\\widehat{\\text{Peer Adoptions}}$',
  'peer.any.sign.100' = '$\\widehat{\\text{\\# Signs}}$',
  'peer.any.vis.100' = '$\\widehat{\\text{\\# Visible}}$'
)

# Preview Table 7 --------------------------------------------------------------
etable(ivsign, ivvis, ivonlysign, ivonlyvis, ivboth,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE,
       dict = varnames, drop = c("peers.100"),
       signif.code = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
       fitstat = c("n", "ivfall", "sargan", "wh"))

# Save Table 7 ------------------------------------------------------
style_tex <- style.tex(tablefoot.value = '',
                       model.title = '',
                       depvar.title = '',
                       var.title = '\\midrule',
                       stats.title = '\\midrule')

etable(ivsign, ivvis, ivonlysign, ivonlyvis, ivboth, tex = TRUE,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE,
       dict = varnames, drop = c("peers.100"),
       signif.code = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
       fitstat = c("n", "ivfall", "sargan", "wh"),
       style.tex = style_tex,
       file = "output/tables/table_7.tex", replace = TRUE)

